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Abstract 

The vahie-at-risk of a delta-gamma approximated derivatives portfoho can be 
computed by numerical integration of the characteristic function. However, 
while the choice of parameters in any numerical integration scheme is paramount, 
in practice it often relies on ad hoc procedures of trial and error. For normal 
and multivariate ^-distributed risk factors, we show how to calculate the neces- 
sary parameters for one particular integration scheme as a function of the data 
(the distribution of risk factors, and delta and gamma) in order to satisfy a 
given error tolerance. This allows for implementation in a fully automated risk 
management system. We also demonstrate in simulations that the method is 
significantly faster than the Monte Carlo method, for a given error tolerance. 

Key words: value-at-risk, delta-gamma approximation, Fourier inversion, 
characteristic function, error bounds 



1. Introduction 

Value-at-risk calculations of large derivatives portfolios are an important 
part of risk management in many financial institutions. The calculations are 
often performed by the Monte Carlo method and can be very time consuming 
— the main bottleneck is the repricing of the portfolio for different simulated 
market scenarios, in particular when the portfolio contains exotic contracts that 
can only be priced with finite difference schemes or even simulations. A common 
way to speed up computations is to consider a delta-gamma approximation of 
the portfolio: the change in the portfolio value over a fixed time horizon is 
approximated by the quadratic function 

AV ^ S'AS + ^AS'TAS, 

where AS is a vector of changes in the underlying risk factors, and d and F 
are the usual first and second order sensitivities. These are typically available 
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"for free" since most often they are already computed internally for hedging 
purposes. Under certain assumptions on the distribution of the risk factors AS, 
the characteristic function 4>{u) = E[e*"'^^] has a closed form expression. One 
can write the cumulative distribution function as 

P(AV^<., = l-£|ge-.<i.. (1) 

and several authors have proposed to compute the value-at-risk by evaluating 
this (or some related) integral numerically, thus avoiding M onte Carlo simula- 
tions in the value-at-risk computations, see iRouvine i (119971 ) and lDuffie fc Pan! 



( 200lh . In this paper, we describe in detail how to perform this computation, 



taking into account that the parameters in the numerical integration scheme 
must be chosen as a function of the data (that is, the distribution of the risk 
factors, and S and F). 

Any numerical integration scheme involves choosing certain parameters and 
for a particular case it is typically not hard to find suitable parameters by 
trial and error. This might not be good enough in an automated system: one 
would like the parameters to be calculated as an internal part of the algorithm, 
preferably with a guaranteed upper bound on the numerical error. 

There is substantial literature on Fo urier series methods fo r co mputing tail 
probabilities, see the extensive reviews by Abate &: WhittI ( 1992f l and Waller et al. 



and the r eferences t herein . Explicit error bounds are however rare: we 



use a result from iHughettI (|l998h . which allows for computing the probability 



in Equation ([T]) with a guaranteed upper bound e on the error, where e > is 
specified by the user. 

To apply the result, one must find certain constants that bound the behavior 
of 4> and the cumulative distribution function. We show how to compute these 
constants as a function of the data when the risk fac tors (i) are norma l ly dis - 



tributed, and (ii) have a multivariate i-distribution. iGlasserman et al.l (j2002r ) 
derive a closed form expression for the characteristic function in the latter, 
heavy-tailed case. We also perform some simulation experiments, where the 
Fourier method is benchmarked against the Monte Carlo method. This is an 
important comparison to make, since the Monte Carlo method is much faster 
for the delta-gamma approximated portfolio than for the original portfolio (it 
avoids the pricing bottleneck), and it is easy to understand and implement. 

The paper is organized as follows. In Section 2, we review the delta-gamma 
approximation and recall the closed-form expressions for the characteristic func- 
tions in the case of normally a nd multivariate t-distributed risk factors. In Sec- 
tion 3, we state the result from lHughett ( 1998f ) on how to perform the numerical 



integration with explicit error bounds, and discuss how to choose the necessary 
parameters. Section 4 contains numerical experiments, and Section 5 concludes. 

2. The delta-gamma approximation of a portfolio 

In this section we define the delta-gamma approximation of a portfolio, and 
give the closed form expressions for the characteristic function under two differ- 
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ent assumptions on the distribution of the risk factors. 

Let S = {Si, . . . , SpY denote a vector of p risk factors to which a portfoho is 
exposed, and let AS* denote the change in 5' from the current time to the end 
of the horizon At. Let V{S, t) denote the value of the portfolio at time t and 
risk factors S. The delta-gamma approximation to the change V{S + AS, At) — 
V{S,0) in the portfolio value is given by 



AV :^ 9At + S'AS + ^AS'TAS*, 



where 



dt' dS^ ''~dS,dSf hJ- I, 



and all partial derivatives are evaluated at {S, 0). The aim is to compute VaR^, 
the level-7 value-at-risk of AV^, under certain assumptions for the distribution of 
AS. The level-7 value-at-risk is defined as the 7-quantile of the distribution of 
AV: VaR-y := P~^{'y), where P{x) := ¥{AV < x) is the cumulative distribution 
function of AV^. Without loss of generality, we look only at the case 9 = 0: if 
P(Ay - eAt < x) = 7, then ¥{AV <x + 9 At) = 7. 

In the case of normally distributed risk factors AS, the characteristic func- 
tion of AV is known on closed form. 

Proposition 1. Assume that AS ^ A/'(0, S) /or some positive definite matrix 
E. Let Ai, A2, . . . , Ap he the eigenvalues ofT,T, and let A be the diagonal matrix 
with these eigenvalues on the diagonal. There is a matrix C satisfying CC = E 
and CTC = A. Let b = C'6. Then the characteristic function corresponding to 
P is given by 

2 P 

(j){u) = E[e™^^] = e"^ ^^=1 ''^'t^tx-^ 11(1 " iX^u)''^/'^, for u G M. 

The moment generating function is given by ipiu) = 4>{—iu) provided uXj < 1 
for j = l,...,p. 

This result is certainly not new, but the proof is constructive and useful for 
implementat ion!!] We give a proof in the Appendix. 

iGlasserman et al.l ( 20021 ) relax the assumption of normal risk factors and look 



at a heavy-tailed distribution for AS". They assume that AS has a multivariate 
^-distribution with degrees-of-freedom parameter v: the density of AS* is given 

by 



^ The distribut i on of quadratic forms of random variables is extensively studied by 
iMathai fc Provost] | |1992D . 
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for some positive definite matrix S. Here, and only here, r(-) denotes the 
gamma- function; elsewhere throughout the paper, T is the matrix of second 
order sensitivities in the delta-gamma approximation of the portfolio. All the 
marginal d istributions of AS are ^ -distributions with degrees-of-freedom pa- 



rameter V. iGlasserman et alj (j2002l) show that, given x G R, P{x) = ¥{AV < 
x) = Fx(0), where is another distribution function, and derive a closed form 
expression for the characteristic function corresponding to . 

Theorem 2 (Glasserman, Heidelberger and Shahabuddin). Denote the 
eigenvalues of EF by Xi, X2, ■ ■ ■ , Xp and let A be the diagonal matrix with these 
eigenvalues on the diagonal. There is a matrix C satisfying CC — E and 
C'TC = A. Let b = C'S. Then P{x) = Fx{Q), where the characteristic function 
corresponding to the distribution Fx is given by 

</.(u) = (l-2e(«))-''/2f[(l-^A,«)-V^ 

with 

1 P 1 
lux 1 Y^,- 



The moment generating function is given by ip{u) — (j){~iu), provided XjU < 1 
for j = 1, . . . ,p and £^{—iu) < i. 



Remark 1. IGlasserman et alJ (|2002f) also consider a generalized copula model, 
where the marginal distributions of AS* are t-distributions with possibly different 
degrees of freedom. The generalization comes down to modifying S and F in a 
straight forward way, and it applies here too. 



3. Calculating VaR-^ by Fourier inversion 

Let F denote the cumulative distribution function for some continuous dis- 
tribution, and let (f) denote the corresponding characteristic function. Assume 
that 4> is known on closed form, but that F is unknown. For a fixed a; S M, 
(|Hughetti il998. Theorem 10) tells us how to compute F{x) with desired accuracy 
by numerical integration of <j). In this section we recall that result and discuss 
how to apply it in the setting of the previous section for computing VaR^. 

Theorem 3 (Hughett). Suppose (i) that there exists constants A and a > 1 
such that F{-y) < and 1 - F{y) < A\y\-" for all y > 0, and (ii) that 

there exist constants B and /3 > such that \4>{u)\ < B\u /2tt\~P for all u G R. 
Then, for constants < I < 2/3, T > and N > 0, the distribution function 
F{x) may be approximated by the truncated Fourier series 

N/2-1 

9ix):^l + 2 Y: Re{G[k]e^'-'^/^), 

k=l 
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where Re(-) denote the real part, and 



1 - cos(27r/fc) 
^= ^^fc H-2nk/T). 

The approximation error on the interval \x\ < lT/2 is bounded by 
\F{x)-g{x)\ < C{f3+l,N/2) + AT-''Li{l,a), 

TT 

where ^ denote the Hurwitz zeta functioi^ and 

L^{1, a) (1/2)-" + 2C{a, 1 - ^0 + C(a, 1 + ^0 + C(a, 1 - ^0- 
Furthermore, for any D > and e > 0, choosing I, T and N such that 

Q<1<\ and rLi{l,a) <2"+\ (2) 
D , „ 2/3A^'/" 



and 



T>- and T>- -) , (3) 



N >2 + 2T 



enf] 



suffices to guarantee that \F{x) — g{x)\ < e for every \x\ < D/2. It is always 
possible to choose I to meet the given conditions. 

Remark 2 (Hughett). The conditions on I, while comphcated, involve only 
the value of a. It is thus feasible to precompute the smallest possible values 
of I, for selected values of a, using a numerical zero-finding algorithm. Table 1 
give the optimal values of I for selected values of a. 

Remark 3. Using Theorem [3J it is straight forward to compute VaR-y both in 
the case of normally and multivariate i-distributed risk factors. For normally 
distributed risk factors, take (j) from Proposition [T] and look for a solution to 
g{x) = 7 with any standard zero-finding method. For multivariate i-distributed 
risk factors we similarly look for a solution to gx{0) = 7, where gx{0) — .g(0) is 
calculated from Theorem [3] with from Proposition [2l where (p depends on x. 
In both cases we end up with x* such that \P(x*) — 7I < e. 

Theorem [3] takes as input parameters A and a, related to the unknown dis- 
tribution function F, and B and /?, related to the known characteristic function 
(j). The rest of this section discusses how to choose these parameters in the 
setting of Section [2l 



■^The Hurwitz zeta function is defi ned by (^{z, a) := "^"^—lik + a) " . For computing it with 
desired accuracy, see l lHuehett]|l998l . Lemma 4). 
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Q 


^opt 


1.125 


0.0855 


1.25 


0.1874 


1.5 


0.3530 


2 


0.4666 


3 


0.4955 


4 


0.4991 


5 


0.4998 


10 


0.5000 



Table 1: Optimal choice of i as a function of a (in Theorem [Sjl, from iHughetti lll998t) . 



3.1. How to choose A and a 

The computational effort in using Theorem [3] is essentially proportional to 
{3A/ey^°', so ideally A and a > 1 should minimize this quantity subject to 

F{-y) < A\y\^" and 1 - F{y) < A\y\-" for aU y > 0, (4) 

where F ^ P for normal risk factors, and F = for the multivariate t- 
distribution. This problem might be difficult to solve, since the constraints 
involve the unknown distribution function F. It is however not important to find 
a truly optimal solution: it suffices to make (SA/e)^/" as small as possible with 
little computational effort. If the moment generating function ip corresponding 
to F is finite in a neighborhood of zeroU one can use Chernoff's bounds, see 
Chernofj (|l952l) : F{~y) < e^"^V(-w) and 1 - F{y) < e""^V(w) for all y and 



all M > 0. A sufficient condition for A and a to satisfy the constraints (|4]) is 
thus that Ay~°' > tp{u)e~'^y , for all y > and some u > 0, where 

V'(u) := max{?/'(M), 

Consider for a moment u > as given, and assume that ili{u) < oo. The 
following proposition gives the solution to the optimization problem 

^ min^^ f — j s.t. Ay-"" > V'(^i)e""^ for all y > 0. (5) 

Proposition 4. The solution to the optimization problem is given by a{u) := 
log(3'0(u)/£) and A(u) :— iJj{u)e~°'^"\a{u)/u)"^"K The minimal function value 
is \og[2>^{u)/e)/u. 

The proof is found in the Appendix. A tractable way of computing a (sub- 
optimal) minimizer of (3A/e)^/" that satisfies the constraints (jl]) is to take 
A = A(u*) and a = a(w*), where u* solves 

minlog(3i/' (6) 

u>0 

This is a one dimensional optimization problem that can be solved with standard 
methods. 



•^If ij) is finite in a neighborhood of 0, then tl){u) = <f)(—iu), so this function is also known 
on closed form. 
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Remark 4. For normal risk factors, the optimization problem ([S]) needs to be 
solved only once. The resulting A and a can be used to compute P{x) for any 
X e R. For the multivariate ^-distribution, P{x) — Fx{0), so the optimiza- 
tion problem needs to be solved once for each x, each time using the moment 
generating function corresponding to F^. 

3.2. How to choose B and (3 

The computational effort in using Theorem [3] is essentially proportional to 
{6B/Trepy^^ , so for fast computations one would like to choose B and /? that 
minimize this quantity, subject to \(j){u)\ < B\u/2n\^^ for all u. In the general 
case, solving this problem comes down to analyzing the characteristic function. 
For the two special cases of normally and multivariate t-distributed risk factors, 
we have the following bounds. The proof is found in the Appendix. 

Proposition 5. Let Ai, . . . , Ap be the eigenvalues of SF, and let I denote any 
subset o/{l,...,p}. Let I3j := |/|/2 and Bj :^ (27r)-^^ Hje/ l-^jT^^^- V ^t^ «s 
given either by Proposition[l\ or Proposition\^ then \(j){u)\ < Bj\u/2n\~'^' . 

A tractable strategy for choosing B and (3 is to solve 

/ \Tref3iJ 

where /?/ and Bj are defined in Proposition [SI and the minimum is taken over 
all the subsets of {1, . . . , p}. 

Proposition 6. Let {ji, . . . be an ordering of {I, . . . ,p}, such that \Xj-^\ > 
|Aj2 1 > • ■ ■ > |Ajp I . Denote Lk := {ji, . . . , jfe} for k = 1, . . . ,p, and let L denote a 
solution to the optimization problem Then I — Ik for some k G {1, . . . 

Remark 5. By Proposition [SI the optimal subset / (and thus /?/ and Bj) can 
be found by computing (6-B/^/7re/3/^)^/'^^fc for k — l,...,p and choosing the 
minimum. Note that this needs to be done only once, even in the case of 
multivariate t-distributed risk factors — although (f> depends on x, the optimal 
parameters Bj and f3i do not. 

4. A large scale numerical experiment 

In this section the value-at-risk for a large, simulated options portfolio is 
calculated with two different methods: (i) the Fourier method from the previous 
section, and (ii) the Monte Carlo method. The results are then compared with 
respect to speed and accuracy. 
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e 



N 



time 



10 
10 
10 
10 



,-6 



216 
314 
442 
774 



0.12 s 
0.15 s 
0.21 s 
0.29 s 



Table 2: Normally distributed risk factors. The integer N from Theorem[3]and the total time 
for the Fourier method computation of VaR-, as a function of e, for 7 = 0.01. 



Simulating a large options portfolio 

We simulate a portfolio that models a real portfolio held by a market maker 
in vanilla options: it contains 10^ vanilla options written on 30 different un- 
der lyings. Each option is long or short with probability ^, and a call or a put 
option with probability ^ . The maturity of each option is drawn from a uniform 
distribution on [10/252, 1] (10 days-1 year), and the moneyness for each option 
is drawn from a truncated normal distribution with mean 1 and standard devia- 
tion 0.1, bounded below by 0.5 and above by 1.5 — the options are between 50% 
out-of-the-money and 50% in-the-money. The returns from the underlyings are 
drawn (i) from a normal distribution or (ii) from a multivariate f-distribution 
with ;/ = 5 degrees of freedom. In both cases, the covariance matrix is AtE, 
where E is a 30 x 30 diagonal matrix with each diagonal element drawn from 
a uniform distribution on [0.1^,0.3^] and At :— 10/252 (so we are looking at 
the 10 day value-at-risk) . The nominal value for each option is drawn from a 
uniform distribution on [10"*, 10^]. 

Having simulated the portfolio, we compute the deltas and gammas for each 
option. Finally, the deltas and gammas are summed up for each underlying, to 
give S and F. 

4-2. Performance of the Fourier method 

Given E, S and F from the simulation, we compute Ai, . . . , A30 and 61, ... , 630 
in Theorems [I] and [21 As seen in the proof of Proposition [1] this is done by 
computing the Cholesky factorization of S, and diagonalizing the matrix SF. 
This gives us the complete closed-form expression for the characteristic function 
and the moment generating function. Let £ > be given, and compute the 
necessary parameters. A, a, B and /3, as described in Remarks|4]and[5l The level- 
7 value-at-risk is then computed as described in Remark [3l The algorithm gives 
a number x* , which is close to VaR^ = P~^(7) in the sense that \P{x*)—^\ < e. 

The computational effort is proportional to the integer N, so we report this 
together with the total computational timeQ for various choices of e in Table 2 
(normal risk factors) and Table 3 (multivariate ^-distributed risk factors). 

4-3. Performance of the Monte Carlo method 

To compute an estimate x* of VaR-y — P^^{'y) by the Monte Carlo method, 
simulate M samples from the A/'(0, S)-distribution (or the multivariate i-distribution). 



*The computations were performed in the software R on a standard workstation. 
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e 




N 


time 


10" 


a 


1000 


0.43 s 


10" 


4 


1508 


0.65 s 


10" 




2190 


0.92 s 


10" 


6 


3070 


1.28 s 



Table 3: Multivariate f-distributed risk factors. The integer N from Theorem [3] and the total 
time for the Fourier method computation of VaR-y as a function of e, for 7 = 0.01. 



compute the corresponding values of AV and order these to get AV(i) < ■ • • < 
AV(Af). Then x* := AV(|-^/^-|), where [M7] is the smallest integer greater 
than or equal to Mj. The computed x* corresponds to some confidence level 
7* :— P{x*), hopefully not too far from 7. To study how large sample size M 
we need to reach a desired level of accuracy we look at the number e such that 
P(|7 — 7*1 < e) — 1 — p for some small positive number p. 

For a moderate sample size M, the corresponding e can be estimated experi- 
mentally in the following way: run the Monte Carlo R times, order the resulting 
estimates x^^^ < ■ ■ • < a;^^^ and compute the corresponding confidence levels 
7*1) < •■• < 7(*it) by evaluating 7?^^ := P{x*^f,s), for A: = 1, . . . , i?, using the 
Fourier method with an extremely low error. Then take 



1 



7(r(l-p/2)fll) -l{[pR/2i) 



where [pi?/2j denotes the largest integer smaller than or equal to pR/2. 

For a large sample size M, we can compute e from the approximate distribu- 
tion of 7*. Since a;* is the [M7]th order statistic, it is asymptotically (as M 
00 ) normally distributee!! with mean VaR-y and variance 7 ( 1 — 7) / (MP' ( VaR.^ ) ^ ) . 
If P is approximately linear in a neighborhood of VaR^, then 7* = P{x*) is 
approximately normally distributed with mean 7 and variance 7(1 — 7)/Af by 
Gauss' approximation formulas. So, s ~ $^-^(1 — p/2)-\/7(l — l)/M, where $(•) 
is the cumulative distribution function of the standard normal distribution. 

Figure [1] show e estimated experimentally for M — 10'^, 10'*, lO'^, with R = 
500 and p = 0.01. It also shows e computed theoretically from the approximate 
distribution of 7*, for M = 10^, lO'', . . . , 10^ — the two methods for computing 
£ agree nicely for M > 10^. If e is small enough, we can thus use the theoretical 
approximate distribution of 7* to compute the sample size M needed to make 
IPdT^T*! < e) ~ 1— p for a given e — see Table 3 for p = 0.01 and various choices 
of £. The total computing time for M = 6.5 • 10** is 1.7 s for normally distributed 
risk factors, and 2.4 s for multivariate t-distributed risk factors. It is moderately 
more time consuming to sim ulate from the mul t ivaria te ^-distribution than from 
the normal distribution, see Glasserman et al. I (|2002l) . 



This is true if Ay has a continuous non-zero density at P"H7)> see [Arnold et al.l lll992l') . 
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log^o(M) 



Figure 1: The graph shows £ as a function of the Monte Carlo sample size M, where £ is 
computed experimentally (stars) and from the approximative theoretical distribution of 7* 
(circles). 



10 



M 



time (normal) time (mult, t) 



10 
10 
10 



1-5 



6.5 • 10* 
6.5 • lO" 
6.5 • 10* 



1.76 s 2.45 s 

3 min 4 min 

4.5 h 6.5 h 



Table 4: The number of Monte Carlo samples M needed to make P(|7* — t| < e) > 0.99 as a 
function of £. We also give the computation times for both the normal and the multivariate 
t-distributions. 



4-4- Comparison of the Fourier method and the Monte Carlo method 

Comparing the results from the Fourier method and the Monte Carlo method, 
we see that the Fourier method is superior both with respect to speed and ac- 
curacy. For £ = 10~^ and normally distributed risk factors, the Fourier method 
takes 0.12 seconds to compute x* such that I7 — P{x*)\ < e, while the Monte 
Carlo method takes 1.76 s to compute x* such that P(|7 - P{x*)\ < e) > 0.99. 
As e ^ 0, the Monte Carlo method becomes impractical: a 10 times smaller 
£ gives 100 times longer computing time. The Fourier method however scales 
nicely: for e ~ 10~^ the whole computation still takes only about 0.3 s. The 
picture is similar for multivariate t-distributed risk factors: the Fourier method 
is significantly faster than the Monte Carlo method, especially for very small e. 

5. Conclusion 

We have given a complete description of how to compute the valuc-at-risk 
of a delta-gamma approximated portfolio by numerical integration of the char- 
acteristic function, in the case of normally or multivariate t-distributed risk 
factors. In particular, we showed how to calculate the necessary parameters in 
the numerical integration scheme as a function of the data (the distribution of 
the risk factors, and S and F), in order to satisfy a given error tolerance. Thus 
we avoid ad hoc choices of parameters, and allow for using the method in a 
fully automated risk management system. Our numerical experiments illustrate 
that the method is significantly faster than the Monte Carlo method — we con- 
clude that the Fourier method is a highly competitive alternative to the Monte 
Carlo method for computing the value-at-risk of delta-gamma approximated 
portfolios. 
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Appendix 

Proof of Proposition[l] Let L be the Cholesky factor of E, so LL' = E and 

L~^AS is a vector of independent, standard normals. Since E is a symmetric 
matrix, L'TL has real eigenvalues, and these eigenvalues are the same as those 
of LL'T — EP, namely Ai,...,Ap. Moreover, L'TL = OAO' , where O is an 
orthogonal matrix whose columns are eigenvectors of L'PL. Since O'L'TLO = A 
and {LO){LOy = LL' — E, setting C = LO produces the required matrix. 

Introduce the identity matrix / LOO' L^^ = {L^^yOO' L' in the expres- 
sion for AV: 

AV = S'AS + ^AS'TAS 

= 5' LOO' L-^ AS + ]^AS'{L-'')'00'L'TLOO'L~^AS 

= {0'L'5)'{0'L-^AS) + ]^{0'L-^AS)'{0'L'TLO){0'L-^AS). 

Write AS = (A5i, . . .,ASp)' := O'L-'^AS and recall that b ^ C'S ^ O'L'S. 
Then 

Ay ^ b'AS + -As'aAS = J2 bj^Sj + 2 H ^i^^J ' 

where ASi,...,ASp are independent, A/'(0, l)-distributed random variables. 
Let J := {j G {1, . . . ,p}; Xj ^ 0} and write 



AV = 
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where the second equahty arises from completing the squares. Define Qo := 
j:^^jbjASj - l/2E,gj&j/A,, and + b^/X,)^, for j G J. Then 

AV — Qo + ^/'^^j^j XjQj is a linear combination of independent random 

variables — Qo and Qj, j G J, are independent since AS*!, . . . , ASp are inde- 
pendent, so the characteristic function of AV is given by 



(f>{u) = E[e™'^^] = E[e™^°] [| E|e 



^^J(Aj/2)«Qj1 

i6J 



The distribution of Qq is normal with mean — 1/2 X]je / ^j/'^^i variance 
a.nd Qj has a non-central x^-distributior0 with 1 degree of freedom 
and non-centrality parameter &^/A^, for j g J. If follows that 




iXjU 



□ 



Proof of Proposition [4l Assume that iIj{u) < oo and consider for a 
moment a > 1 as given. Then the constraint will be satisfied if A is taken large 
enough. For the smallest possible A, the graph of the function f{y) — Ay~°' 
will not intersect the graph of h{y) = tlj{u)e^'^y , but be tangent to it at some 
point y* . This is expressed in the following system, 

f(.y*) = h{y*), 
fiy*) = h'iy*), 
f"{y*) > h"{y*), 

which is solved by y* = a/u and A{a^u) = %l){u)e^°'{a/u)'^ . Plugging the 
expression for A(a, u) into the minimization problem ^ reduces it to 
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a I 3-0(it) \ 



This target function has the unique minimum a{u) := log{3tp{u)/e). The result 
follows, since (3A(a(u), = i log(3V;(u)/£). 



□ 



^The characteristic function of a x^-distributed random variable Z with k degrees of free- 
dom and non-centrality parameter a is given by E[e'"^] = e" (i _ i2u)~''/'^. 
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Proof of Proposition [5j If </) is given by Proposition [U then 
Mn)Mn), with Mn) exp (-i ^^^^ ^^(^) j^p^^^^ 

is given by Proposition [21 then (f>{u) — 4>^{u)4>2{u) with (jy^iu) : 
(1 - 2$(u))-''/2. Since ^ iAju) = y+I^ + *TTI^' foUows that 
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e J < 1. 



It also follows that 



|l"2C(u)| = 



1 ^ 



+ A2u2 



+ iIni(-2^(M)) 



> 1, 



so |<?!)3(m)| = |1 - 2^(u)|-'-'/2 < 1. 

To bound |02(u)| > AjM|, first note that |1— mAj| = y'l + A^m^ > niax{l, |Aju|}, 
so |1 - iAjii|^^/2 < niin{l, \\ju\^^/'^}. This implies that 



-1/2 



i=i 



for any subset I c{l,...,p]. It follows that |(/)('u)| < Sl'u/27r|-'3 with (3 = \I\/2 
and B = {2tt)-P H^g/ lA^ T^/^ both for (j> = ^!)i(/>2 and </> = (/>302- 

□ 

Proof of Proposition [6l Assume the contrary to the Proposition: I ^ Ik 
for k = 1, . . . ,p. Then there exists m, m' e {1, • • ■ to < to', such that 
jm ^ I and jrn' e I. Let /' (/ \ {j™'}) U {jm}- Then /?/ = /?// and, unless 
lAj™ I = |Aj„,, I, < Bj — this contradicts the optimality of /. The degenerate 
case |Aj„J — only arises from the trivial ambiguity in ordering numbers 

of the same size. To rule it out, we adopt the following convention: if s indices 



jt+i,- ■ ■ ,jt+s correspond to |A 



= ••• = |A 



3t+s 



and if / contains fc < s of 



these indices, then it contains the k first ones: jt+i, ■ ■ ■ ,jt+k- 



□ 
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